home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / IGAM.C < prev    next >
C/C++ Source or Header  |  1995-10-03  |  4KB  |  217 lines
  1. /*                            igam.c
  2.  *
  3.  *    Incomplete gamma integral
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double a, x, y, igam();
  10.  *
  11.  * y = igam( a, x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * The function is defined by
  18.  *
  19.  *                           x
  20.  *                            -
  21.  *                   1       | |  -t  a-1
  22.  *  igam(a,x)  =   -----     |   e   t   dt.
  23.  *                  -      | |
  24.  *                 | (a)    -
  25.  *                           0
  26.  *
  27.  *
  28.  * In this implementation both arguments must be positive.
  29.  * The integral is evaluated by either a power series or
  30.  * continued fraction expansion, depending on the relative
  31.  * values of a and x.
  32.  *
  33.  *
  34.  *
  35.  * ACCURACY:
  36.  *
  37.  *                      Relative error:
  38.  * arithmetic   domain     # trials      peak         rms
  39.  *    DEC       0,30         4000       4.4e-15     6.3e-16
  40.  *    IEEE      0,30        10000       3.6e-14     5.1e-15
  41.  *
  42.  */
  43. /*                            igamc()
  44.  *
  45.  *    Complemented incomplete gamma integral
  46.  *
  47.  *
  48.  *
  49.  * SYNOPSIS:
  50.  *
  51.  * double a, x, y, igamc();
  52.  *
  53.  * y = igamc( a, x );
  54.  *
  55.  *
  56.  *
  57.  * DESCRIPTION:
  58.  *
  59.  * The function is defined by
  60.  *
  61.  *
  62.  *  igamc(a,x)   =   1 - igam(a,x)
  63.  *
  64.  *                            inf.
  65.  *                              -
  66.  *                     1       | |  -t  a-1
  67.  *               =   -----     |   e   t   dt.
  68.  *                    -      | |
  69.  *                   | (a)    -
  70.  *                             x
  71.  *
  72.  *
  73.  * In this implementation both arguments must be positive.
  74.  * The integral is evaluated by either a power series or
  75.  * continued fraction expansion, depending on the relative
  76.  * values of a and x.
  77.  *
  78.  *
  79.  *
  80.  * ACCURACY:
  81.  *
  82.  *                      Relative error:
  83.  * arithmetic   domain     # trials      peak         rms
  84.  *    DEC       0,30         2000       2.7e-15     4.0e-16
  85.  *    IEEE      0,30        60000       1.4e-12     6.3e-15
  86.  *
  87.  */
  88.  
  89. /*
  90. Cephes Math Library Release 2.0:  April, 1987
  91. Copyright 1985, 1987 by Stephen L. Moshier
  92. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  93. */
  94.  
  95. #include "mconf.h"
  96.  
  97. #define BIG  1.44115188075855872E+17
  98. extern double MACHEP, MAXLOG;
  99.  
  100. double igamc( a, x )
  101. double a, x;
  102. {
  103. double ans, c, yc, ax, y, z;
  104. double pk, pkm1, pkm2, qk, qkm1, qkm2;
  105. double r, t;
  106. double lgam(), exp(), log(), fabs();
  107. double igam();
  108. static double big = BIG;
  109.  
  110. if( (x <= 0) || ( a <= 0) )
  111.     return( 1.0 );
  112.  
  113. if( (x < 1.0) || (x < a) )
  114.     return( 1.0 - igam(a,x) );
  115.  
  116. ax = a * log(x) - x - lgam(a);
  117. if( ax < -MAXLOG )
  118.     {
  119.     mtherr( "igamc", UNDERFLOW );
  120.     return( 0.0 );
  121.     }
  122. ax = exp(ax);
  123.  
  124. /* continued fraction */
  125. y = 1.0 - a;
  126. z = x + y + 1.0;
  127. c = 0.0;
  128. pkm2 = 1.0;
  129. qkm2 = x;
  130. pkm1 = x + 1.0;
  131. qkm1 = z * x;
  132. ans = pkm1/qkm1;
  133.  
  134. do
  135.     {
  136.     c += 1.0;
  137.     y += 1.0;
  138.     z += 2.0;
  139.     yc = y * c;
  140.     pk = pkm1 * z  -  pkm2 * yc;
  141.     qk = qkm1 * z  -  qkm2 * yc;
  142.     if( qk != 0 )
  143.         {
  144.         r = pk/qk;
  145.         t = fabs( (ans - r)/r );
  146.         ans = r;
  147.         }
  148.     else
  149.         t = 1.0;
  150.     pkm2 = pkm1;
  151.     pkm1 = pk;
  152.     qkm2 = qkm1;
  153.     qkm1 = qk;
  154.     if( fabs(pk) > big )
  155.         {
  156.         pkm2 /= big;
  157.         pkm1 /= big;
  158.         qkm2 /= big;
  159.         qkm1 /= big;
  160.         }
  161.     }
  162. while( t > MACHEP );
  163.  
  164. return( ans * ax );
  165. }
  166.  
  167.  
  168.  
  169. /* left tail of incomplete gamma function:
  170.  *
  171.  *          inf.      k
  172.  *   a  -x   -       x
  173.  *  x  e     >   ----------
  174.  *           -     -
  175.  *          k=0   | (a+k+1)
  176.  *
  177.  */
  178.  
  179. double igam( a, x )
  180. double a, x;
  181. {
  182. double ans, ax, c, r;
  183. double lgam(), exp(), log();
  184. double igamc();
  185. static double big = BIG;
  186.  
  187. if( (x <= 0) || ( a <= 0) )
  188.     return( 0.0 );
  189.  
  190. if( (x > 1.0) && (x > a ) )
  191.     return( 1.0 - igamc(a,x) );
  192.  
  193. /* Compute  x**a * exp(-x) / gamma(a)  */
  194. ax = a * log(x) - x - lgam(a);
  195. if( ax < -MAXLOG )
  196.     {
  197.     mtherr( "igam", UNDERFLOW );
  198.     return( 0.0 );
  199.     }
  200. ax = exp(ax);
  201.  
  202. /* power series */
  203. r = a;
  204. c = 1.0;
  205. ans = 1.0;
  206.  
  207. do
  208.     {
  209.     r += 1.0;
  210.     c *= x/r;
  211.     ans += c;
  212.     }
  213. while( c/ans > MACHEP );
  214.  
  215. return( ans * ax/a );
  216. }
  217.